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ABSTRACT 

Context. The fidelity of Smoothed Particle Hydrodynamics (SPH) simulations of accretion disks depends on how 
Artificial Viscosity (AV) is formulated. 

Aims. We investigate whether standard methodology is reliable in this regard. 

Methods. We test the operation of two methods for selective application of AV in SPH simulations of Keplerian Accretion 
Disks, using a ring spreading test to quantify effective viscosity, and a correlation coefficient technique to measure the 
formation of unwanted prograde alignments of particles. 

Results. Neither the Balsara Switch (B) nor Time Dependent Viscosity (TDV) work effectively, as they leave AV active 
in areas of smooth shearing flow, and do not eliminate the accumulation of alignments of particles in the prograde 
direction. The effect of both switches is periodic, the periodicity dependent on radius and unaffected by the density of 
particles. We demonstrate that a very simple algorithm activates AV only when truly convergent flow is detected and 
reduces the unwanted formation of prograde alignments. The new switch works by testing whether all the neighbours 
of a particle are in Keplerian orbit around the same point, rather than calculating the divergence of the velocity field, 
which is very strongly affected by Poisson noise in the positions of the SPH particles. 
Conclusions. 
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1. Introduction 

In Smoothed Particle Hydrodynamics (SPH) simulations 
which involve strongly convergent flows, it is essential to in- 
clude Artificial Viscosity (AV) (Monaghan & Gingold 1983; 
Monaghan & Lattanzio 1985; Monaghan 1992). It causes 
close, approaching SPH particles to repel each other, with 
a force which increases with approach velocity. The result 
is that particles in colliding streams are rapidly decelerated 
and the streams do not pass through each other. Such par- 
ticle interpenetration would be an unphysical situation and 
cannot be allowed in a simulation. It is important to note 
that AV was never intended to replicate the behaviour of 
real viscosity (Monaghan 2005). 

SPH is now used widely to simulate the evolution of 
disks (Bate et al 2003; Lodato & Rice 2004; Rice et al 2003a, 
2003b, 2004; Mayer et al 2004; Schaefer et al 2004; Boffin 
et al 1998; Watkins et al 1998a and 1998b, Stamatellos 
et al 2008 and 2009). The behaviour of viscosity in real 
circular shear flow in fluids is understood analytically and 
experimentally (eg Feynman 1964), and the evolution of 
accretion disks is predicted to be dominated by effects such 
as turbulence which behave like viscosity, causing angular 
momentum to be transferred outwards, enabling matter to 
fall in and accrete onto a central object (Shakura & Sunyaev 
1973; Lynden Bell & Pringle 1974). It is often assumed that 
turbulence is generated and moderated by the magneto- 
rotational instability (Balbus & Hawley 2002). 

However, AV, if allowed to operate unchecked in a sim- 
ulated Keplerian Disk, does not operate only when con- 
vergence occurs, but also acts to slow down neighbouring 
particles as they overtake because of the differential rota- 



tion. The magnitude of AV is therefore, systematically, too 
high. Clarke (2009) has shown that it is a key requirement 
for correct modelling of the evolution of disks that numer- 
ical viscosity is kept very low, otherwise it dominates the 
pseudo viscous forces, such as self-gravity, in the simula- 
tion. It is therefore crucial to disable AV in normal, equi- 
librium Keplerian shear flow, and to activate it only when 
it is needed to capture a shock. 

In this paper we set out to test the Balsara Switch 
(B), and Time Dependent Viscosity (TDV), for their effec- 
tiveness in disabling AV in non-convergent Keplerian shear 
flow. In Sect. 2 we explain the theoretical operation of B 
and in Sect. 3 we examine the actual values of B calcu- 
lated for a snapshot of a smoothly rotating accretion disk. 
In Sect. 4 we investigate the variation of B with time for in- 
dividual SPH particles. In Sect. 5 we turn to TDV and look 
at actual values calculated for a simulated Keplerian disk 
in which there is no convergence. In Sect. 6 we look at the 
spreading of a ring of particles, as a quantitative measure of 
the effectiveness of B and TDV in reducing viscous shear. 
We also quantify, for the first time, the previously reported 
tendency for AV to cause SPH particles to form prograde 
alignments. A new approach to controlling AV in disks, by 
simply checking whether neighbours are in Keplerian orbit 
around the same object, is explained and evaluated in Sect. 
7 and these results are discussed in Sect. 8. 



2. The Balsara Switch 

The Balsara Switch (B) is a well known tool for enabling 
the rapid application of AV in the presence of a shock 
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(Balsara 1989), and works as follows. A multiplicative 
'Balsara factor', 



Balsara Values 



Bi = 



V-v 



V x v 



(1) 



is calculated for each SPH particle i, where v is the velocity 
field and the derivatives are evaluated at position r^. 

The viscous acceleration between each pair of par- 
ticles is then multiplied by the mean of the values of B, 
thus 



(Bi 



Bj) 



n, 



(2) 



In a Keplerian disk, with all matter orbiting in cen- 
trifugal balance about a mass M, V x v | and V • v | at 
radius r are given by : 



(GM\ 1/2 

Vxv = [ipr) 

I V-v I = 0. 



(3) 



| V x v | is therefore finite while | V ■ v | is zero, and hence 
B is zero for all r. However, if particles start to converge at 
a particular location, V • v | increases, so B increases, ap- 
proaching a limiting value of 1.0, when | V-v |^| V xv |. 
Multiplying the viscosity coefficients by B should effectively 
switch off the AV when a disk is rotating stably with a 
Keplerian velocity profile, but switch it on, selectively, for 
any particles in an area of true convergence. 

3. Noise problems with the SPH implementation 

of B 

In SPH simulations B is calculated for each particle using 
Eqn. [TJ substituting numerical estimates V • v |sph and 
V x v I sph for the true values of V ■ v and V x v : 



V-v. 



i |SPH 



V X V, I SPH 




(4) 
(5) 



From Eqn. [31 the value of B should be zero for a stably 
rotating disk of particles in Keplerian orbits, where there 
is no convergence. However, as Fig. [T] demonstrates, this 
does not happen in practice in SPH simulations. 

Figure [1] shows that high values of B were calculated 
for 80,000 SPH particles in purely Keplerian orbits around 
a star mass = IMq. The two dimensional ring was 
initialised with randomly positioned particles, such that 
the surface density was Gaussian around a maximum value 
at radius = 10AU. An area of convergence was created 
around coordinates (7,7), by adding convergent velocities 
up to 1.5km/s to the Keplerian velocities, which are around 
150km/s at these radii. In Fig. [TJa quadrant of the ring was 
viewed face-on, the particles orbiting in the anti-clockwise 
direction. Particles with B > 0.3 were plotted as green dots, 
and those with B > 0.6 as red dots. Clearly B was found to 
be high, not only in the converging region, but all over the 
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Fig. 1. Values of B for a portion of a ring of particles in 
Keplerian orbits. 

disk, with some values as high as 0.8. The particles with 
high B were not scattered randomly over the disk, rather 
they appeared to form alignments, predominantly in the 
'leading' direction. 

The neighbours of a particle in a Keplerian disk fall 
into four quadrants of approaching or receding particles. 
Particles ahead and in outer orbits and particles behind 
and on inner orbits are approaching, so Vy ■ is nega- 
tive. Particles ahead on inner orbits and behind on outer 
orbits are receding, so Vy- • r.y is positive. The calculation 
of V • v | sph adds together positive values of v^-iy, from 
receding neighbours and negative values from approaching 
ones. If neighbours are evenly distributed, V • v |sph is 
small, as the positive and negative components cancel out. 
However if neighbours are unevenly distributed, Poisson 
noise causes V • v |sph to increase in magnitude. 

This is shown in column 2 of Table []] where the mean 
value of B was calculated for two disks. The first disk con- 
tained randomly positioned particles in perfect Keplerian 
orbits, and B should therefore be zero. An average value 
of B — 0.11 was consistently found in this and all other 
randomly populated disks, for numbers of particles rang- 
ing from 80000 to 1 million. The second, settled, disk was 
created by allowing the first to evolve using pressure and 
AV, controlled by B. There was no convergence in this disk, 
and because the particles were more evenly spaced, B was 
closer to zero. B fell to less than half the value of the orig- 
inal disk, to B = 0.052, and remained steady. Column 3 of 
the same table shows the mean value of B calculated for a 
convergence zone created within the same settled disk. 



4. Variation of V • v |sph with time 

In a Keplerian disk, all particles are constantly overtak- 
ing their outer neighbours and being overtaken by inner 
neighbours. As a result, three or more particles often find 
themselves in radial alignment. In this case, these aligned 
neighbours are in opposite quadrants, so they are all re- 
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Table 1. Mean value of B and MK, for 80,000 SPH particles 
orbiting in a Keplerian disk. 



neighbour does, which takes a time 
the distance travelled is therefore 



-. During this time 



Disk B B MK NIC 

Status Disk Conv Disk Conv 

Random 0.11 0.003 

Settled 0.052 0.116 0.060 0.774 




0.5 




Time yr 



Fig. 2. B plotted against time, for SPH particles at three 
different radii in a Keplerian disk around a IMq star. 



ceding or all approaching. They therefore contribute to a 
large negative value of V • v |sph while they are approach- 
ing radial alignment, and then a high positive value, once 
they have passed radial alignment and are spreading out. 
This explains the orientation of the linear features observed 
in Fig. [TJ Chance alignments of particles will produce high 
values of B for any orientation of alignments, if B is calcu- 
lated for both positive and negative values of V • v |sph- 
However if B is set to zero for positive values of V • v |sph 7 
high B values only occur for points in alignments with outer 
particles ahead and inner particles behind. 

Figure [5] shows the value of B for three SPH particles 
in different orbits within a Keplerian disk. As B has been 
set to zero when V • v |sph is positive, we see the positive 
half of three periodical functions. The solid line indicates a 
particle at radius 1.4AU, the dashed line a particle at 3 AU 
and the dotted line a particle at 10AU. Orbital periods for 
these radii are 1.6y, 5y and 30y. The discontinuities in B 
are due to modifications to the neighbour lists, which affect 
the calculated values of V • v |sph and V x v |sph- Clearly 
B is varying quasi-periodically, with a frequency decreasing 
with the radius of the orbit. 

This is explained by considering the 'overtaking' speed 
of a neighbour. Consider a particle orbiting mass M at 
radius r, and its neighbour orbiting at radius r + h. The 
overtaking orbital speed, v re i, of the inner particle is 



Vrel 



dr 



(GM) 



1/2 



2r 3/2 



(6) 



In order to overtake the neighbour, we allow the par- 
ticle to travel a distance h further along its orbit than its 



• v orb = 2r 

Vrel 



(7) 



which is approximately one third of an orbit. 

If the inner particle travels 2r along its orbit, the outer 
must travel 2r + 2h, just to keep its station. So, adding 
this effect to that calculated in Eqn. [?1 we have the particle 
moving 3h past a neighbour a distance h away in the time 
it takes to travel one third of the way round its orbit. This 
should be approximately the period of one cycle of the value 
of B. 

The variation of V-v |sph 7 and hence of B, with time is 
therefore predicted to be periodic with frequency approxi- 
mately three times the orbital frequency. A particle orbiting 
at radius 1.4 AU with an orbital period of 1.6 year should 
show 4 cycles of B in a 2 year period. For 3 AU the period 
is about 5 years and we expect just over one cycle. For 10 
AU, with an orbital period 30 years, we only expect to see 
about g of the cycle. This agrees well with what is seen in 

Fig. m 

Once an aligned structure has formed, all the particles 
in that structure will have high V • v |sph, values, and 
therefore both B and the magnitude of the AV (which 
is also calculated using V • v |sph), will be large for all 
these aligned particles. They will therefore tend to stick 
to their neighbours more than other particles. The longer 
such structures exist, the more angular momentum they 
will transfer. Imaeda and Inutsuka (2002) observed that 
linear structures are a pronounced feature of SPH simula- 
tions of shear flow which use AV. This raises the worrying 
possibility that AV itself could be enhancing the formation 
of filaments in disk simulations, with B doing nothing to 
moderate this tendency. This could lead to misleading re- 
sults. 



5. Time Dependent Viscosity 

A second method for selectively switching on AV in shocks 
while keeping it switched off at other times, was proposed 
by Morris and Monaghan (1997), and is referred to as Time 
Dependent Viscosity (TDV). 

Each SPH particle i is given its own viscosity parameter 
CKi, which evolves with time, according to the equation 



dai 
~~dt 



where 5, is the source term, 



Si = max(-V • v |sph 



,0). 



Tj is the c-folding time, which is set to be 



CiCi ' 



(8) 



(9) 



(10) 



where hi is the smoothing length for viscous and pressure 
forces, Cj the sound speed and C\ a parameter chosen to 
control the rate at which oti decays back to a* . 

Clearly a, increases if V-v |sph is negative, and decays 
to a minimum value a* otherwise. This provides a rapidly 
increasing viscosity around a particle which is entering a 
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shock, where V • v |sph is large and negative for a number 
of successive timesteps. Once the shock is passed, V-v |sph 
falls to small or even positive values, and the viscosity dies 
away again. Morris and Monaghan (1997) suggest that set- 
ting C\ to 0.2 should ensure that the viscosity dies down 
to the minimum value within about 5 smoothing lengths. 

This method is not helpful for preventing spuriously 
high viscosity in a Keplerian accretion disk, the problem 
being that once again V • v |sph is being used as an indi- 
cator of convergence. As we have shown, V • v |sph varies 
periodically with a frequency dependent on radius, as par- 
ticles overtake one another, causing alignments to form. It 
is impossible to tune the TDV method to filter out such low 
frequency noise without making it much too slow to react 
to shocks. 



6. AV Tests 

6.1. Spreading ring test 

As a simple test of the magnitude of effective shear vis- 
cosity in a real SPH simulation, with and without B or 
TDV, a ring spreading test was carried out. Lynden-Bell 
and Pringlc (1974) predicted that a ring of particles will 
spread under the influence of viscosity, the rate of spread- 
ing increasing with viscosity. Murray (1996) and Flebbe 
et al.(1994) conducted SPH simulations to demonstrate 
this effect, starting with a ring initialised with a Gaussian 
Surface Density profile. 

The SEREN SPH code (Hubber 2010, in prep) was 
used. This is a code specifically developed for the study of 
star formation. A ring of 80,000 SPH particles, total mass 
O.IMq was created, with a Gaussian distribution of surface 
density about a maximum at radius 10 AU. This was set 
in Keplerian orbit about a 1M Q star, with no self gravity 
within the ring and no pressure forces, and allowed to evolve 
only under the influence of the star's gravity and AV, for 
a fixed simulated real time duration (200 years). Multiple 
Particle Time Stepping was used with a maximum range of 
10 4 from the shortest to the longest timesteps. As the ring 
evolved, all the particles were binned in 100 bins by ra- 
dius, the mean surface density calculated for each bin, and 
the maximum value stored as T, max (t). AV was calculated 
for 50 neighbours, and only applied for negative values of 
V • v | sph- The viscosity parameter a was varied from 0.1 
to 1.0, while j3 was kept at 0.0. This was then repeated with 
B implemented, as defined in Eqn.[T]and again with TDV. 
The results are shown in Table [H Column 1 gives the value 
of AV a parameter, column 2 the method used to moder- 
ate AV, ( B, TDV, A/7C or NONE). Column 3 gives the 
measured reduction in maximum surface density, <5£(%). 

The rings did not evolve smoothly, the central density 
oscillating as the ring evolved. However, there was a mea- 
surable trend, and the first four rows of Table [2] show that 
the ring spread more quickly, with the central density falling 
more rapidly, with increasing values of a. 

With no convergence within the ring, dE should be 0.0 if 
B and TDV were effective in switching off AV in shear flow. 
Table [2] shows that implementing B and keeping a = 1.0 
did decrease the effective viscosity, giving a smaller <5X than 
that obtained when a was reduced to 0.1. By comparison, 
our mean B value for 50 neighbours in Tabled] was 0.052. 
TDV worked slightly less well, being larger than for B. 



Table 2. Reduction in maximum surface density, <5E(%) found 
for various AV prescriptions in the ring spreading test. 



a 


AV Switch 


<5E(%) 


0.0 


None 


0.00 


0.1 


None 


1.08 


0.5 


None 


2.18 


1.0 


None 


3.39 


1.0 


BS 


0.72 


1.0 


TDV 


0.92 


1.0 


NK. 


0.00 



Points in alignment with close neighbours 




x AU 

Fig. 3. Prograde (red crosses) and retrograde (black cir- 
cles) particle alignment in a disk evolved with pressure and 
viscous forces. 

6.2. Measurement of particle alignment 

We have suggested that prograde alignments of particles 
will persist because they are held together by AV, which is 
artificially high in these alignments. The data from the ring 
spreading experiments was therefore analysed to quantify 
any such effect. Each particle in turn was examined, to- 
gether with any close neighbours, defined as lying within 
0.5h, where h is the SPH smoothing length. The coordi- 
nates of the close neighbours were projected into radial 
and tangential components, within the disk, and if at least 
three points were found within 0.5/i the correlation coef- 
ficient cc was calculated for the group of points. A value 
of cc > 0.5 indicated a linear alignment in the prograde 
direction, while cc < —0.5 indicated a linear alignment in 
the retrograde direction. Figure [3] illustrates that this cal- 
culation correctly identified the prograde and retrograde 
alignments (red crosses and black circles respectively) in a 
portion of the disk, rotating anticlockwise. 

Prograde and retrograde alignments were then counted, 
and the results are shown in column 7 of table [31 In the 
disk illustrated in Fig. [3l 24% of particles were in prograde 
alignments, and 10% in retrograde, a prograde surplus of 
14%. 
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Table 3. Percentages of prograde (8A + ) and retrograde (5 A ) 
alignments found in disks evolved with various prescriptions of 
AV. 
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None 


IS 
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n i 


No 
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None 


Z8 
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1Z 


l.U 


INO 


None 


Zo 


11 


1 A 
l.U 


JNO 


h 


z ( 


1 /I 

14 


1.0 


No 


TDV 


28 


13 


1.0 


No 


NIC 


19 


18 


0.0 


Yes 


None 


20 


20 


0.1 


Yes 


None 


24 


15 


0.5 


Yes 


None 


22 


15 


1.0 


Yes 


None 


35 


11 


1.0 


Yes 


B 


30 


12 


1.0 


Yes 


TDV 


32 


12 


1.0 


Yes 


NIC 


23 


17 



Initial, random, particle arrangements had 18% parti- 
cles in prograde and 18% in retrograde alignments, so there 
was no excess in prograde alignments in the initial disk. The 
first four rows of Table[3]show that there were also no excess 
prograde alignments when the disk was evolved with no AV 
(a = 0.0), but when AV was applied, the number of pro- 
grade alignments increased to 26 — 28% while the number 
of retrograde alignments decreased to 11 — 14%. The next 
three rows show that using B resulted in 27%/14% pro- 
grade excess, and the TDV AV 28%/13%. Thus, although 
both methods moderated the spreading effect of AV, the 
build up of an excess of prograde alignments was still large. 

The trial was repeated with the addition of pressure 
forces. This increased the number and imbalance of align- 
ments. With a = 0.0, no effective AV, the number of align- 
ments was seen to increase in both directions to 20%. As 
AV was increased the number of prograde alignments in- 
creased and retrograde decreased, the largest imbalance 
being seen for a = 1.0, when 35% of particles were in pro- 
grade alignments, and only 11% in retrograde. B and TDV 
produced only slightly fewer alignments than this. Thus 
pressure forces were seen to exacerbate the problem of pro- 
grade alignments being formed preferentially while retro- 
grade alignments either did not form or were disrupted. 

7. Keplerian recognition approach 

Another way of ruling out areas of apparent convergence 
which arc in fact simply regions of Keplerian shear, is 
to compare neighbours' positions and velocities with a 
Keplerian template. Here we use the facts that in a 
Keplerian disk, the velocity is always perpendicular to the 
position vector of a particle relative to the centre of rota- 
tion, and the quantity v 2 r is constant over all of the disk. 
Clearly, the position of the centre of rotation must be known 
to apply these tests. 

For each particle pi, we first calculated v\ri for the parti- 
cle itself. We then examined each of the Nn ei b neighbours 
of the particle. A neighbour was counted as anomalous if 
its velocity was not perpendicular to its position vector, 
indicated by rj.Vj/\rjVj\ > e. We also counted the point 



as anomalous if vfri/v^rj > 1.0 + e or < 1.0 — e. The er- 
ror value e yielded good discrimination when set to a value 
of 0.03. Large e failed to identify important discrepancies 
from Keplerian flow, while very small e caused viscosity to 
be switched on when pressure forces caused slightly non- 
Keplerian velocity gradient. 

The number of anomalous particles Na were accumu- 
lated and then divided by the total number of neighbours 
to give a viscosity switch NIC = Na/Nj^sib- NIC was 
therefore near zero for regions where all particles were in 
Keplerian orbit around the same point, and increased to 1 
as more and more neighbours fail to conform to a Keplerian 
pattern. Table [T] shows that in a disk of 80000 particles or- 
biting with a perfectly Keplerian velocity profile, the mean 
value oi NIC was 0.003. When the disk was allowed to set- 
tle, using both pressure and AV, NIC increased to 0.06. 
Although this mean value is slightly greater than the mean 
value of B, the distribution of high values was very differ- 
ent. High values of B were registered all over the disk, while 
high values of NK were mainly found in the edge regions. 

NIC was found to be effective at activating AV in areas 
of convergence. This is illustrated in Fig. HJ where NIC 
is compared with B. A ring was evolved with particles in 
Keplerian orbits with pressure and AV activated, and then 
a circular region of converging velocities was added to the 
Keplerian flow. All fifty neighbours of a point were given 
a convergence velocity proportional to their distance from 
the central point and pointing towards it, with the max- 
imum convergence velocity being 1% of the orbital veloc- 
ity of the central particle. This was added to the orbital 
velocity of each particle. Note that this was quite a stiff 
test, the magnitudes of the artificial converging velocities 
being of the same order of magnitude as those due to the 
Keplerian velocity shear within the convergence zone, and 
much smaller than the sound speed. The particles in the 
region of the convergence are shown, with crosses to indi- 
cate particles for which B or the non-Keplerian indicator 
A/7C, were greater then 0.2, green indicates values higher 
than 0.5, blue greater than 0.8. It can be seen that both 
methods correctly indicate the area of convergence, but 
NIC increases for all particles in the convergence region. 
B is raised for most particles in the convergence zone, but 
is also incorrectly raised for alignments of particles which 
are not in the convergence zone. The mean value of NIC 
within the convergence zone is 0.78, and .06 outside while 
mean B is only 0.116 inside the zone and 0.052 outside. 

The ring spreading test was repeated using the NIC fac- 
tor, and the results are presented in Table [5J NIC kept AV 
switched off and the ring did not spread at all. There was 
also only a very small increase in the growth of prograde 
linear alignments (from 18 to 19%) when the NIC factor was 
used. When pressure forces were added to the simulation, 
the NIC factor did not completely suppress the formation of 
alignments, but performed much better than B and TDV. 

Finally, to test the NIC method with a more realis- 
tic simulation, an isothermal disk was evolved contain- 
ing a Jupiter-mass object. Initial conditions were isother- 
mal, T=10K, surface density £ = Z(lAU)(R/au)- J / 4 , 
M disk = O.O1M . M star = 1M Q , The disk extended from 
5 to 15 au and planet mass was 1 Mj at radius 10 au. 
500,000 SPH particles were used and the simulation ran for 
470yrs, about 15 orbital periods of the planet. The results 
are shown in Fig. [5j It can be seen that NIC was effective 
at capturing the shocks, and the predicted spiral wave was 
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Fig. 4. B (top) and MK. (bottom) values for points in a 
Keplerian disk, with an area of converging flow superim- 
posed at coordinates (-5.6,7.6). 



observed. For comparison, the same simulation was per- 
formed using B, and the results are similar, but the spirals 
and shocks are slightly better defined when using the NIC 
method. 



8. Discussion 

The SPH estimate of the divergence of viscosity | V • v | , 
(Eqn. 01 when used in a Keplerian Disk, is dominated by 
the Poisson noise arising from the imperfectly spaced SPH 
particles (Cartwright et al 2009). As the estimate of | V ■ v | 
is the key factor in calculating both B and TDV, both 
of these methods erroneously register the apparent conver- 



gence of overtaking neighbours as a signal that AV should 
be applied. 

B is prone to an alignment artifact, varying periodically 
with a frequency ~ 3 times the orbital frequency, which is 
the overtaking frequency for neighbours. This is indepen- 
dent of h and therefore independent of the number of SPH 
particles used. 

TDV is also not ideal for use in Keplerian disks if used 
with V ■ v | sph, as a source term. V • v |sph is the source 
of the low frequency alignment artifact which compromises 
B, and applying a low pass filter to it does not remove the 
effect. 

AV if applied, even at quite a low level, in a Keplerian 
Disk, results in the formation of alignments of particles 
in the prograde direction. This can be quantified by cal- 
culating the correlation coefficient cc of close neighbours 
of particles, positive values cc > 0, 5 indicating prograde 
alignments, cc < —0.5 indicating retrograde alignments. A 
randomly populated disk of particles has no excess of ei- 
ther, having about 18% of each. The balance is unaffected 
if the disk is evolved with only pressure forces and no AV. 
However, as soon as AV is used, even at low levels using 
B, the particles start to form an excess of 12 — 16% pro- 
grade alignments. Using pressure forces in addition to AV 
causes an even larger surplus of alignments to appear, up to 
24% being measured when maximum AV was implemented. 
This should be of great concern for the disk-modelling com- 
munity. Normally, disk simulations result in the formation 
of spiral arms in the retrograde direction, which may then 
form gravitational instabilities. The action of AV to form 
and preserve prograde alignments, while destroying retro- 
grade ones, may well be interfering in this process. 

In order to prevent particle inter-penetration, AV 
must be applied promptly when areas of convergence are 
detected. B is not as effective here as might be ex- 
pected. When an area of convergence is superimposed upon 
Keplerian shear, the mean B value within the shear zone is 
only B = 0.116, whereas values greater than this are seen 
elsewhere in the undisturbed regions of the disk. 

A solution to the problem of keeping AV switched off in 
regions of Keplerian Shear, and reaching the very low lev- 
els pointed out as necessary by Clarke (2009), may lie in a 
pattern recognition approach. Here we demonstrate a very 
simple method, which checks that all neighbours of a point 
are travelling orthogonally to the radius of the disk, and 
that the value of v 2 r is constant within a small tolerance 
for all neighbours. This simple method results in AV being 
switched off consistently in areas of smooth Keplerian shear 
(mean AfJC — 0.06), while switching on more effectively in 
areas of simulated convergence (AftC = 0.77 in convergence 
zone). It also markedly reduces the formation of prograde 
alignments of SPH particles. In a disk evolution simulation 
with pressure and AV enabled, B and TDV resulted in 22% 
and 24% surplus prograde alignments, while the excess with 
the JV7C switch was only 6%. The A/7C method requires the 
orientation and centre of rotation of the disk to be known, 
but this is true for most disk simulations, and easily de- 
duced for others. The calculation is no more complicated 
than that for the B, so there is no increased computational 
overhead, and it is very amenable to parallelisation. A/7C 
also has the interesting feature of being variable, by alter- 
ing the value of e. A high value of e ~ 0.1 will keep AV 
switched off in normal thermal motions superimposed on 
Keplerian shear, but will only switch on for particle i if 



A. Cartwright: A switch for Viscosity in SPH simulations 



7 




J 



Fig. 5. Disks incorporating a planetary mass object, evolved using left, Balsara Switch and right the AftC factor. 



neighbouring particles j have non-Keplerian approach ve- 10. Acknowledgements 



locities Vij greater than 10% of Vi. If e ~ 0.01, AV will 
be applied to smaller discrepancies from Keplerian motion, 
which may include some thermal motions within the disk. 

Note that MK, does not identify convergence, merely 
the absence of Keplerian flow. However, as it is then used 
in conjunction with AV, which does detect convergence, 
the two methods in combination are effective in applying a 
decelerating viscous force only when SPH particles are in 
truly convergent flow within a Keplerian disk. 

Finally, it should be noted that the flaw in the calcu- 
lation of | V • v | also applies in the calculation of | V • a | . 
Large divergence in the acceleration or velocity fields can 
apparently be detected when in fact the variation is due to 
noise in the particle positioning. Both quantities should be 
used with caution, particularly as diagnostic indicators in 
the creation of sink particles. 

9. Conclusion 

Both B (Balsara 1989) and the TDV approach (Morris and 
Monaghan 1997) to controlling AV in circular shear flow 
fail for the same reasons: SPH estimates of V • v appear to 
detect convergence in steady shear flow. AV controlled by 
these methods is therefore applied in all areas of a simu- 
lated Keplerian accretion disk, with results which will be 
unprepresentative of real viscous evolution, and compro- 
mise the results (Clarke, 2009). AV results in the selective 
formation of prograde alignments of SPH particles, which 
can be measured using a correlation coefficient method. We 
find that a disk populated with randomly placed SPH par- 
ticles initially has 18% particles in prograde alignments and 
18% retrograde. After being allowed to evolve using pres- 
sure forces and AV, these figures change to 28% and 12%. 
Clearly the use of AV changes the arrangement of parti- 
cles in the experiment, even when very low levels of AV are 
applied. 

An alternative method of controlling AV, based on iden- 
tifying Keplerian velocity characteristics, shows promise, 
being more effective both at switching on in regions of con- 
vergence and off in pure Keplerian flow, and significantly 
reducing the accumulation of excess prograde alignments of 
particles. 
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